Removal of alpha-wave artefacts in MEG data by independent component analysis 
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1 Introduction 

It is well known from electro-encephalographic 
(EEG) and magneto-encephalographic (MEG) 
recordings that the human brain generates spon¬ 
taneous oscillatory signals at various frequencies. 
Probably the best known are the a waves with a 
typical frequency of 10 Hz, which are related to 
reduced vigilance or a state of low mental activity. 
These a waves can severely corrupt the stimulus 
triggered average in evoked response experiments. 
Motivated by the observation of 10 Hz contamina¬ 
tion in averaged Visually Evoked Fields (aVEF) ob¬ 
tained without inter-stimulus interval randomization, 
we performed an Indepedent Component Analysis 
(ICA) on raw unaveraged data. It was the aim to 
identify the source of the 10 Hz contamination and 
to remove it. The Temporal Difference source SEPa- 
ration (TDSEP) algorithm [1,2] was applied as it is 
very well suited for oscillatory data containing sig¬ 
nals with distinct spectral features [3]. The identifica¬ 
tion of noise sources is generally an important topic 
as their efficient removal will improve the analysis of 
individual evoked epochs and their variability. 

Other work [4] has studied the relation between a 
waves and motor responses using a different ICA al¬ 
gorithm. 

2 Methods 

2.1 Visually evoked fields 

In normal subjects the MEG was recorded with a 
planar 49 channel gradiometer [5] centered over T6. 
Grey tone images of houses and faces were presented 
for 500 ms and masked for another 1.5 s leading to 
a total epoch length of 2.0 s. Using a block design 
either a face or a house was the attended target (press 
button) to study attention related effects. Here we 
focus on the separation of artefactual signal compo¬ 
nents. 

The four field maps shown in Fig. 1 correspond to 
the most important features of an aVEF calculated 
from 90 non-consecutive epochs. Image presentation 
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Figure 1: Four field maps of a typical aVEF measured 
over T6. The image presentation started at 0 ms. The 
field pattern changes rapidly as the map sequence 
shows. 


started at 0 ms and a first maximum is reached 160 ms 
later as shown in the first map. The data in Fig. 1 cor¬ 
respond to the case that the attention was focused on 
faces and a face was presented. 

In Fig. 2a) we show the superimposed time traces of 
the 49 MEG channels for the same aVEF as in Fig. 1. 
The main activity starts at 160 ms. Weaker activity 
is already present at 130 ms and even at -100 ms, be¬ 
fore image presentation at 0 ms. This activity cannot 
be related to the image presentation without violating 
causality. Spectral analysis of the aVEF reveals that 
the signal before image presentation is mainly an os¬ 
cillation with a frequency of 10 Hz. A spectral anal¬ 
ysis of the raw unaveraged data gives a strong 10 Hz 
peak, which increases in strength as the experiment 
proceeds. The aVEF shown in Fig. 2a) is an aver¬ 
age from the last block of an experiment containing 
in total three blocks each of 20 minutes duration. We 
have chosen the last block as the 10 Hz oscillations 
are strongest in this block and therefore the need to re¬ 
move the contamination is most urgent. Certain sub¬ 
jects showed strong 10 Hz oscillations throughout the 
experiment. 



a) 
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Figure 2: a) Butterfly plot of the aVEF shown in 
Fig. 1. Even before image presentation at 0 s large 
oscillations with a periodicity of 100 ms are visible, 
b) From the data in a) the component ICA1 of Fig. 3 
has been removed, c) The components ICA1, 3, 5, 10 
and 11 (not shown) have been removed from the data 
in a). 

2.2 Temporal Difference source SEParation 

In MEG recordings an unknown number of brain and 
external (noise) sources is present. It is difficult to 
model this in full generality and commonly a simple 
linear model is used. It is assumed that the recorded 
data are a linear superposition of n sources s % ( t ) with 
a constant source pattern aji, which are measured in a 


sensor system having m channels yielding Xj ( t ) time 
traces. To describe this in an invertible linear equation 
m > n is required and the model can be written as 

x(f)=As(f). (1) 

The sources s, (f) and their patterns are generally un¬ 
known. A demixing matrix W is searched for, which 
fulfills 

u(f)=Wx(f), (2) 

where the uflt) are identical to the sflt) disregarding 
an indeterminancy in scale and ordering. This needs 
to be resolved in the sense, that the sources can be or¬ 
dered by their strength to discriminate between strong 
and weak sources. This is achieved by normalizing 
the vectors of W and then calculating the Euclidean 
norm (L2-norm) of the resulting uflt). 

ICA is a method to estimate the demixing operator W 
using the assumption of statistical independence be¬ 
tween sources. This independence is equivalent to a 
probability density function factorising into the indi¬ 
vidual densities as shown: 

= IlPi(si)- (3) 

The factorization means, that the covariance matrix 
and time-delayed covariance matrices are diagonal 
with respect to the sflt). The TDSEP [1, 2] algorithm 
utilizes this as it tries to simultaneously diagonalize 
several time-delayed covariance matrices to estimate 
the demixing matrix W. 

For interpretation and assessment the blindly calcu¬ 
lated Ui(t) and their patterns W -1 can be compared 
afterwards to secondary information. A simple ex¬ 
ample is the heartbeat. Together with the MEG a sin¬ 
gle electrocardiogram (ECG) channel can be recorded 
and this time series can be compared to each of the 
uft). 

3 Results 

3.1 TDSEP applied to raw data 

The TDSEP algorithm was applied to a dataset of 80 s 
length corresponding to 40000 datapoints for each of 
the 49 sensor channels. Twenty time delayed covari¬ 
ance matrices offset by 2 ms were used. This dataset 
is 10 % of the full data set. The results do not dif¬ 
fer considerably using the full set. By using smaller 
data subsets it is possible to compare the result from 







which are severely contaminated by eye blinks. 

The result of the decomposition using TDSEP is 
shown in Fig. 3. For each component three types of 
visualization are given, the field map (i.e. the columns 
of W _1 ), the time series (i.e. ufit)), and the spec¬ 
trum of the time series. The time series is only a short 
excerpt from the 80 s used in the calculation. The 
components are ordered by the value of their F2-norm 
given as ’F2’ in the figure. 

The first component ICA1 shows amplitude modu¬ 
lated 10 Hz oscillations. In ICA2 the lower frequency 
components are nearly as strong as the 10 Hz oscil¬ 
lations. The field maps of ICA1 and 2 resemble the 
aVEFs in Fig. 1 at 280 and 240 ms. The smooth and 
almost monopolar field map of ICA4 implies that it is 
related to cortical activity, although it cannot be iden¬ 
tified with parts of the aVEF. Component ICA4 is ei¬ 
ther due to a failure of the ICA model (e.g. n > m, 
statistical dependence between the s t (t), etc.) or due 
to the presence of a real physiological signal not visi¬ 
ble in the average and not known from secondary in¬ 
formation. 

Component ICA 5 is related to the QRS-complex of 
the heartbeat as its time series is similar to a simul¬ 
taneously recorded ECG. The raw data averaged with 
respect to the R-peak showed for certain times a field 
map very similar to ICA 5 and was for other times 
similar to ICA 3. The field maps of 7CA 3 and 5 show 
gradients typical for a remote source as it is the heart. 
The distortions in the gradient are due to the detec¬ 
tor geometry relevant for remote sources. Component 
ICA10 has a field maximum at the lower edge of the 
sensor. This was the position of the neck in the mea¬ 
surement and muscular tension seems to be a possi¬ 
ble origin. This is supported by its spectrum having a 
broad maximum at 40 Hz. A component due to 50 Hz 
power line noise looks very different (not shown). 


Figure 3: ICA decomposition of multichannel raw 
MEG data using TDSEP. Component ICA1 has a peak 
at 10 Hz in its spectrum and is attributed to a waves, 
ICA2 has a broad spectrum and its map is similar to 
the map at 240 ms in Fig. 1. ICA4 cannot be asso¬ 
ciated immediately to a known source. ICA3 and 5 
are attributed to the heartbeat, and finally ICA10 is 
possibly caused by muscular activity in the neck. 


different regions in the dataset and to exclude regions 
of strong external noise or rest periods of the subject, 


3.2 Noise supression by demixing 

As the TDSEP decomposition in Fig. 3 has resulted in 
two components containing 10 Hz oscillations it ap¬ 
peared promising to demix the measured signal and 
recombine it omitting one or several of the compo¬ 
nents. This is done mathematically by setting the cor¬ 
responding column(s) in W 1 to zero. 

Firstly the component with the strongest 10 Hz signal 
is omitted ( ICA1 ), and the result is shown in Fig. 2b). 
The oscillations before image presentation at 0 ms 
and between 300 and 400 ms are clearly reduced as 
is positive activity at 130 ms and a negative peak at 






255 ms. The physiologically most interesting early 
activity between 160 and 240 ms remains essentially 
unchanged. The reduction of the 10 Hz oscillations 
in the aVEF due to omitting component ICA1 implies 
that the 10 Hz oscillations present in the raw data ap¬ 
pear in the average. 

Omitting 7CA2 the aVEF is strongly modified in the 
region between 160 and 240 ms indicating that 7CA2 
is related to the early aVEF. The dual role of ICA1 and 
2 in describing parts of the aVEF together with oscil¬ 
latory activity is either due to an incomplete demixing 
or it is of physiological origin, i.e. VEF and the 10 Hz 
oscillations originate from the same cortical location. 
In Fig. 2c) the components ICA1, 3, 5, 10 and 77 
(50 Hz power line noise, not shown) have been omit¬ 
ted. The signal to noise ratio before the physiologi¬ 
cally expected aVEF onset at 100 ms seems further 
improved implying artefacts in the aVEF due to the 
heartbeat or power line noise. This has to be assessed 
yet quantitatively. 

The two TDSEP components containing 10 Hz oscil¬ 
lations are the strongest signals of the raw data sec¬ 
tion used in Fig. 3 as they have the largest values of 
F2. Therefore it should be possible to observe them 
in bandpass filtered raw data. Using a filter of width 
5 -15 Hz field maps similar to ICA1 and 2 can be seen 
at the same times as in Fig. 3. Taken all facts together 
we identify the 10 Hz oscillations with a waves. This 
is further supported by the amplitude modulation oc¬ 
curring in the time series of ICA1 and 2. 

The amplitude of the a waves in the filtered data is up 
to 1 pT. The 10 Hz oscillations before stimulus onset 
peak at 50 fT in Fig 2a). Assuming the a waves to be 
uncorrelated to the stimulus they should be reduced 
in the aVEF by a factor of 90 1 / 2 « 9.5. Given the 
intermittent nature of a waves this agrees well with 
observation. 

Note that simple bandpass filtering cannot yield the 
same result as the TDSEP decomposition. The filter 
result is due to the combined pass band spectral den¬ 
sity of all signals, i.e. a waves, heartbeat, VEF and 
other sources. 

4 Conclusions 

Motivated by a 10 Hz contamination in an aVEF we 
performed an TDSEP-ICA on MEG raw data. Sig¬ 
nificant components can be attributed to a waves, the 
heartbeat, muscular tension in the neck and the VEF 
itself. Omitting the strongest component attributed to 
a waves the aVEF is considerably improved. Clearly 


the effectiveness of the method depends on the orien¬ 
tation of the dominant a field map compared to phys¬ 
iologically interesting parts of the aVEF, which was 
very favourable here and for other subjects. Only in 
the case of a subject reporting fatigue after the experi¬ 
ment it was not possible to improve the aVEF. Future 
work will focus on using TDSEP to isolate compo¬ 
nents solely attributable to the VEF. 
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